#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#include <vector>
#include <string>
#include <algorithm>
#include <sstream>
#include <ostream>
#include <iomanip>
using std::vector;
using std::string;
using std::max;

#include <ppl.h>
#include <amp.h>
//#include <amp_math.h>
using namespace concurrency;
//using namespace concurrency::fast_math;

#define FLOAT float

#include "Option.h"
#include "Model.h"
#include "Stopwatch.h"
#include "MTGP-simple/MtgpSimple.h"
#include "Reduce.h"
#include "amp-rng/amp_rand.h"
#include "MTGP-parallel/mtgp32-amp.h"

struct PartialResults
{
    int count;
    double sum, sum2;
    bool timing;
    //Random rng = new Random(0);
    Stopwatch s;

    PartialResults(bool showTiming = true)
    {
        count=0;
        sum=0;
        sum2=0;
        timing = showTiming;
    }
    void add(double oneResult)
    {
        count++;
        sum += oneResult;
        sum2 += oneResult * oneResult;
    }
    string totals()
    {
        double mean = sum / count;
        double var = sum2 / count - mean * mean;
        std::ostringstream str;
        str << std::setw(4)
            << "Trials:" << count 
            << "\tMean: " << (float)mean
            //<< "\tStdDev: " << (float)sqrt(var)
            << "\tStdErr: " << (float)sqrt(var / count);
        if(timing)
        {
            str << "\ttime/ms: " << s.Milliseconds()
                << "\tMtrials/s: " << count/s.Milliseconds()/1000
                << "\t";
        }
        return str.str();
                        
    }
};

//        void CalcMonteCarloSingleThreaded(float rate, float time, European option, int trials)
//        {
//            std::cout << "CalcMonteCarloSingleThreaded" << std::endl;
//            //Random rng = new Random(0);
//            float drift = (rate - 0.5f * option.vol() * option.vol()) * time;
//            float df = (float)exp(-rate * time);
//            int nextOutput = 2;
//            PartialResults results;
//            for (int i = 0; i < trials; ++i)
//            {
//                float U = (float)(rng.Next(1, 100000 - 1) / (float)(100000));
//                float W = NormalCDFInverse(U);
//                float S = option.spot * (float)exp(drift + W * option.vol * (float)Math.Sqrt(time));
//                float C = df * option.Payoff(S);
//                results.add(C);
//                if ((int)Math.Log10(i) > nextOutput)
//                {
//                    nextOutput = (int)Math.Log10(i+1);
//                    Console.WriteLine(results.totals() + option.ClosedForm(rate, time));
//                }
//            }
//        }
//
//        void CalcMonteCarloParallel1(float rate, float time, European option, int trials)
//        {
//            Console.WriteLine("CalcMonteCarloParallel1");
//            //object guard = new object();
//            float drift = (rate - 0.5f * option.vol * option.vol) * time;
//            float df = (float)exp(-rate * time);
//            int nextOutput = 0;
//            PartialResults fullResults;
//
//            Concurrency::parallel_for(
//                0, trials, 
//                () => new PartialResults(),
//                (i, loopState, partialResults) =>
//                {
//                    float U = (float)(partialResults.rng.Next(1, 100000 - 1) / (float)(100000));
//                    float W = NormalCDFInverse(U);
//                    float S = option.spot * (float)exp(drift + W * option.vol * (float)Math.Sqrt(time));
//                    float C = df * option.Payoff(S);
//                    partialResults.add(C);
//                    return partialResults;
//                },
//                (localPartialResults) =>
//                {
//                    lock (guard)
//                    {
//                        fullResults.sum += localPartialResults.sum;
//                        fullResults.sum2 += localPartialResults.sum2;
//                        fullResults.count += localPartialResults.count;
//                        if ((int)Math.Log10(fullResults.count) > nextOutput)
//                        {
//                            nextOutput = (int)Math.Log10(fullResults.count);
//                            Console.WriteLine(fullResults.totals() + option.ClosedForm(rate, time));
//                        }
//                    }
//                }
//            );
//        }
//
//        void CalcMonteCarloParallel2(float rate, float time, European option, int trials)
//        {
//            Console.WriteLine("CalcMonteCarloParallel2");
//            object guard = new object();
//            float drift = (rate - 0.5f * option.vol * option.vol) * time;
//            float df = (float)exp(-rate * time);
//            int nextOutput = 0;
//            PartialResults fullResults;
//
//            var partitioner = System.Collections.Concurrent.Partitioner.Create(0, trials);
//            System.Threading.Tasks.Parallel.ForEach(partitioner,
//                () => new PartialResults(),
//                (range, loopState, partialResults) =>
//                {
//                    for (int i = range.Item1; i < range.Item2; i++)
//                    {
//                        float U = (float)(partialResults.rng.Next(1, 100000 - 1) / (float)(100000));
//                        float W = NormalCDFInverse(U);
//                        float S = option.spot * (float)exp(drift + W * option.vol * (float)Math.Sqrt(time));
//                        float C = df * option.Payoff(S);
//                        partialResults.add(C);
//                    }
//                    return partialResults;
//                },
//                (localPartialResults) =>
//                {
//                    lock (guard)
//                    {
//                        fullResults.sum += localPartialResults.sum;
//                        fullResults.sum2 += localPartialResults.sum2;
//                        fullResults.count += localPartialResults.count;
//                        if ((int)Math.Log10(fullResults.count) > nextOutput)
//                        {
//                            nextOutput = (int)Math.Log10(fullResults.count);
//                            Console.WriteLine(fullResults.totals() + option.ClosedForm(rate, time));
//                        }
//                    }
//                }
//            );
//        }

template<class Model, class Option>
void TestMtgpSimple(MtgpSimple& pmt, const Model& model, const Option& option)
{
    //SETUP
    const extent<2> randomsExtent = pmt._randomsExtent;
    const bool gaussian = true;
    const bool dumpRNG = false;
    const size_t DISPLAYPATHS = 0;

    //DEVICE CODE
    PartialResults results;
    std::cout << "Device:" << std::endl;
    array<float, 2> d_random(randomsExtent);
    pmt.generate_rand_on_amp(d_random, gaussian);
    pmt._acc.wait();
    std::cout << "RNG: " << results.s.Milliseconds() << "ms, Mrng/s=" << d_random.extent.size()/results.s.Milliseconds()/1000 << std::endl;
    if(dumpRNG)
    {
        PartialResults RNGresults1; size_t out=10, c=0;
        for(int n=0; n<d_random.extent[0]; n++) {
        for(int m=0; m<d_random.extent[1]; m++) {
            RNGresults1.add(d_random[n][m]); 
            //std::cout << n << "\t" << m << "\t" << d_random[n][m] << std::endl;
            if(++c==out) {
                std::cout << "RNG " << RNGresults1.totals() << std::endl;
                out*=3;
            }
        }}
    }

    array<float,2> d_path(randomsExtent);
    array<float,2> d_pv(randomsExtent);
    concurrency::parallel_for_each(
        d_random.extent, 
        [=, &d_random, &d_path, &d_pv] 
        (index<2> idx)
        restrict(amp)
    {
        d_path[idx] = model.generatePath(d_random[idx], option.maturity());
        d_pv[idx] = option.payoff(d_path[idx], model.df(option.maturity()));
    });
    double d_sum, d_sumSqr;
    static const int tilrandomsExtent = 16;
    reduction_cascade<tilrandomsExtent, 12800000/16>(
        d_pv.view_as(extent<1>(d_pv.extent.size())), d_sum, d_sumSqr);
    results.count = d_random.extent.size();
    results.sum = d_sum;
    results.sum2 = d_sumSqr;

    std::cout 
        //<< "n=" << pv.size() << "\t"
        << results.totals()
        << "\tBS=" << option.closedForm(model.spot(), model.rate(), model.vol())
        << std::endl;

    std::vector<float> random(randomsExtent.size()), path(randomsExtent.size()), pv(randomsExtent.size());
    copy(d_random, random.begin());
    copy(d_path, path.begin());
    copy(d_pv, pv.begin());
    for(size_t n=0; n<DISPLAYPATHS; n++) 
        std::cout << "n=" << n << "\tr=" << random[n] << "\tp=" << path[n] << "\tpv=" << pv[n] << std::endl;

    //CPU SINGLE CODE
    PartialResults results2;
    std::cout << "CPU single:" << std::endl;
    std::vector<float> s_random(randomsExtent.size());
    pmt.generate_rand_on_cpu(s_random, gaussian);

    if(dumpRNG)
    {
        PartialResults RNGresults; size_t out=10, c=0;
        for(size_t n=0; n<s_random.size(); n++) {
            RNGresults.add(s_random[n]); 
            std::cout << n << "\t" << s_random[n] << std::endl;
            if(++c==out) {
                std::cout << "RNG " << RNGresults.totals() << std::endl;
                out*=3;
            }
        }
    }

    std::vector<float> s_path(randomsExtent.size());
    std::vector<float> s_pv(randomsExtent.size());
    for(size_t idx=0; idx<s_pv.size(); ++idx)
    {
        s_path[idx] = model.generatePath(s_random[idx], option.maturity());
        s_pv[idx] = option.payoff(s_path[idx], model.df(option.maturity()));
    }
    for(size_t n=0; n<DISPLAYPATHS; n++) 
        std::cout << "n=" << n << "\tr=" << s_random[n] << "\tp=" << s_path[n] << "\tpv=" << s_pv[n] << std::endl;

    for(size_t n=0; n<s_pv.size(); n++) 
        results2.add(s_pv[n]);

    std::cout 
        //<< "n=" << pv.size() << "\t"
        << results2.totals()
        << "\tBS=" << option.closedForm(model.spot(), model.rate(), model.vol())
        << std::endl;

    //CPU PARALLEL
    PartialResults results3;
    std::cout << "CPU multi:" << std::endl;
    std::vector<float> h_random(randomsExtent.size());
    pmt.generate_rand_on_cpu(h_random, gaussian);
    std::vector<float> h_path(randomsExtent.size());
    std::vector<float> h_pv(randomsExtent.size());
    concurrency::parallel_for(
        size_t(0), h_pv.size(), size_t(1),
        [=, &h_random, &h_path, &h_pv, &option] 
        (size_t idx)
    {
        h_path[idx] = model.generatePath(h_random[idx], option.maturity());
        h_pv[idx] = option.payoff(h_path[idx], model.df(option.maturity()));
    });
    for(size_t n=0; n<DISPLAYPATHS; n++) 
        std::cout << "n=" << n << "\tx=" << h_pv[n] << std::endl;

    for(size_t n=0; n<h_pv.size(); n++) 
        results3.add(h_pv[n]);

    std::cout 
        //<< "n=" << pv.size() << "\t"
        << results3.totals()
        << "\tBS=" << option.closedForm(model.spot(), model.rate(), model.vol())
        << std::endl;
}

void TestAmpRng(accelerator_view acc, extent<2> randomsExtent)
{
    std::cout << "AMP-RNG" << std::endl;
    extent<1> extent0(randomsExtent[0]);
    tinymt_collection<1> rngs(extent0, acc);
    array<float, 2> d_randoms(randomsExtent, acc);
    Stopwatch s;

    //generate
    parallel_for_each(acc, extent0, [=, &d_randoms] (index<1> idx) restrict(amp)
    {
        auto rng = rngs[idx];
        rng.initialize(123245);
        for(int i=0; i<randomsExtent[1]; ++i)
            d_randoms(idx[0], i) = rng.next_single();
    });
    acc.wait();
    double ms = s.Milliseconds();
    
    //copy back to host
    vector<float> randoms(randomsExtent[0]*randomsExtent[1]);
    copy(d_randoms, randoms.begin());
    double ms2 = s.Milliseconds();

    std::cout << randomsExtent.size() << "(" << randomsExtent[0] << "x" << randomsExtent[1] << ") in " 
        << ms << "ms, Mrng/s = " 
        << randomsExtent.size()/ms/1000 
        << ",\twith copy Mrng/s=" 
        << randomsExtent.size()/ms2/1000 
        << std::endl;
    //for(int i=0; i<5; ++i) std::cout << randoms[i*10000] << " ";
    //std::cout << std::endl;
}

template<class Model, class Option>
void TestMtgpParallel(accelerator_view acc, int size, const Model& model, const Option& option) 
{
    const int loops = 10;
    const size_t DISPLAYPATHS = 1;

    //AMP CODE
    std::cout << "AMP" << std::endl;
    mtgp_amp mtgp(size);
    acc.wait();
    array<float> d_random(size);
    array<float> d_path(size);
    array<float> d_pv(size);

    for(int i=0; i<loops; ++i)
    {
        Stopwatch s;
        mtgp.make_float_random_amp(d_random);
        acc.wait();
        double ms = s.Milliseconds();
        std::cout << size << " in " << ms << "ms, Mrng/s = " << size/ms/1000 << std::endl;

        concurrency::parallel_for_each(
            d_random.extent, 
            [=, &d_random, &d_path, &d_pv] 
            (index<1> idx) restrict(amp)
        {
            d_path[idx] = model.generatePath(d_random[idx], option.maturity());
            d_pv[idx] = option.payoff(d_path[idx], model.df(option.maturity()));
        });
        acc.wait();
        ms = s.Milliseconds();
        std::cout << size << " in " << ms << "ms, MPath/s = " << size/ms/1000 << std::endl;
        double d_sum, d_sumSqr;
        static const int reduction_tile_size = 1024;
        reduction_cascade<reduction_tile_size>(d_pv, d_sum, d_sumSqr);
        ms = s.Milliseconds();
        std::cout << size << " in " << ms << "ms, MResults/s = " << size/ms/1000 << std::endl;
        PartialResults results;
        results.count = size;
        results.sum = d_sum;
        results.sum2 = d_sumSqr;

        std::cout 
            //<< "n=" << pv.size() << "\t"
            << results.totals()
            << "\tBS=" << option.closedForm(model.spot(), model.rate(), model.vol())
            << std::endl;

        std::vector<float> random(size), path(size), pv(size);
        copy(d_random, random.begin());
        copy(d_path, path.begin());
        copy(d_pv, pv.begin());
        for(size_t n=0; n<DISPLAYPATHS; n++) 
            std::cout << "n=" << n << "\tr=" << random[n] << "\tp=" << path[n] << "\tpv=" << pv[n] << std::endl;
    }

    //HOST CODE
    std::cout << "CPU" << std::endl;
    mtgp_amp h_mtgp(size);
    vector<float> h_random(size);
    for(int i=0; i<loops; ++i)
    {
        Stopwatch s;
        h_mtgp.make_float_random_cpu_parallel(h_random);
        //h_mtgp.make_float_random_cpu_singlethread(h_random);
        double ms = s.Milliseconds();
        std::cout << size << " in " << ms << "ms, Mrng/s = " << size/ms/1000 << std::endl;

        vector<float> h_path(size);
        vector<float> h_pv(size);
        concurrency::parallel_for(
            size_t(0), h_pv.size(), size_t(1),
            [=, &h_random, &h_path, &h_pv, &option] 
        (size_t idx)
        {
            h_path[idx] = model.generatePath(h_random[idx], option.maturity());
            h_pv[idx] = option.payoff(h_path[idx], model.df(option.maturity()));
        });
        ms = s.Milliseconds();
        std::cout << size << " in " << ms << "ms, MPath/s = " << size/ms/1000 << std::endl;
        for(size_t n=0; n<DISPLAYPATHS; n++) 
            std::cout << "n=" << n << "\tr=" << h_random[n] << "\tp=" << h_path[n] << "\tpv=" << h_pv[n] << std::endl;

        PartialResults results;
        for(size_t n=0; n<h_pv.size(); n++) 
            results.add(h_pv[n]);
        ms = s.Milliseconds();
        std::cout << size << " in " << ms << "ms, MResults/s = " << size/ms/1000 << std::endl;

        std::cout 
            //<< "n=" << pv.size() << "\t"
            << results.totals()
            << "\tBS=" << option.closedForm(model.spot(), model.rate(), model.vol())
            << std::endl;
    }
    std::cout << std::endl;

}


void main(int argc, char **argv)
{
    int size=256*3 * 100 * 100;
    //auto accelerators = accelerator::get_all();
    //for_each(begin(accelerators), end(accelerators),[=](accelerator acc)
    accelerator acc;
    { 
        std::wcout << std::endl << L"Using device : " << acc.get_description() << std::endl;
        if(!acc.is_emulated)
        {
            //SETUP
            MtgpSimple pmt(extent<2>(128, size/128), acc.get_default_view());
            Forward forward(2.0f, 80.0f);
            European european(2.0f, 90.0f);
            European european2(2.0f, 110.0f);
            BlackScholes blackScholes(100, 0.2f, 0.0f);

            try
            {
                //std::cout << std::setw(4);
                //std::cout << std::endl << "Forward " << forward.maturity() << "\\" << forward.strike() << std::endl;
                //TestMtgpSimple(pmt, blackScholes, forward);
                //TestMtgpSimple(pmt, blackScholes, forward);
                //std::cout << std::endl << "European " << european.maturity() << "\\" << european.strike() << std::endl;
                //TestMtgpSimple(pmt, blackScholes, european);
                //std::cout << std::endl << "European " << european2.maturity() << "\\" << european2.strike() << std::endl;
                ////TestMtgpSimple(pmt, blackScholes, european2);

                for(int i=4; i<65536; i*=4)
                    TestAmpRng(acc.default_view, extent<2>(i, size/i));

                TestMtgpParallel(acc.default_view, size, blackScholes, forward);
            }
            catch(const std::exception& ex)
            {
                std::cout << "Exception: " << ex.what() << std::endl;
            }
        }
    }
    //);
}
